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Current phylogenetic comparative methods generally employ the Ornstein-Uhlenbeck(OU) process 
for modeling trait evolution. Being able of tracking the optimum of a trait within a group of related 
species, the OU process provides information about the stabilizing selection where the population 
mean adopts a particular trait value. The optima of a trait may follow certain stochastic dynamics 
along the evolutionary history. In this paper, we extend the current framework by adopting a rate 
of evolution which behave according to pertinent stochastic dynamics. The novel model is applied 
to analyze about 225 datasets collected from the existing literature. Results validate that the new 
framework provides a better fit for the majority of these datasets. 
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Introduction 


In evolutionary biology, phylogenetic comparative methods (PCMs) are commonly ap¬ 
plied to analyze trait data for a group of species. Since species share evolutionary history, 
a good estimate of a phylogenetic tree, which represents the evolutionary relationship, 
is incorporated in data analysis. In past decades many comparative methods have been 
developed under different evolutionary hypotheses. For instance, the trait of a group of 
related species may rely on a continuous process ranging from either a single Brownian 
motion [^, or Brownian motions (BMv) to an Ornstein-Uhlenbeck (OU) process 
or multiple optima, multiple rate of evolution and multiple strength of selection 
OU process(OUmva) 14. For a more comprehensive review of comparative methods, the 
reader may refer to |24l |. 

Recently some PCMs have been developed for an advanced study of adaptive evolu¬ 
tion in a randomly evolving environment [IS], ll7| . In contrast to the model of correlated 
evolution which sorely predicts the response trait, the model of adaptive evolution esti¬ 
mates the optimal relationship between two traits. To describe the model for adaptive 
evolution, we start with describing the evolutionary behavior for the response trait. 

Let ut be the trait value of a species at time t, Ot be the optimum(evolutionary central 
tendency) of yt, be the rate of adaptation representing the speed of the trait tracking 
on its optimum, (fi be the rate of evolution of yt, and be a white noise. A trait 
value of a species, yt, is a solution of the following Ornstein-Uhlenbeck (OU) stochastic 


* Corresponding author. Email: dcjhwueng@fcu.edu.tw 


1 




August 14, 2015 


JhwuengMaroulas'JAL'submit 


differential equation (SDE) 


dyt = al{6t - yt)dt + (T^dW^. 


( 1 ) 


Several models have been developed and applied widely to analyze trait data adopting 
further assumptions. For instance, the work in 1^ adopts equation ([1]) by assuming 


that the parameters are all time invariant that is a\ = Oy, 6^ = Oy and = ay are 


all constant in equation ([T]). The study in [g] assumes that multiple optima of yt occur 
during the evolutionary process. In this case, 6^ = 0^ where 9^ is a piecewise constant 
value on the time interval where 7 = 1 , 2 , • • • , m and to = 0 is the initial time 

of evolution and tm = T is the time length from t = 0. Their model is then applied to 
study the evolution of the body size of anolis lizard of the northern Lesser Antillean 
where each of these small islands supports either one or two species of anoles of different 
size (one species is large while the other is small). Beaulieu et al. [3] extend the models 
considered in such that the force parameters, at = a^, and rate parameters, at = a^, 
are also allowed to take piecewise constant values on and the models are applied 

to study the genome size evolution within a fairly large flowering plant clade. 

However, the optimum, 6t, typically is not a constant (or piecewise) static parameter. 
Instead its behavior typically evolves according to another independent OU process given 
below, 


de'i = -a<ie\dt + aUw!, ( 2 ) 

where the parameters af,af are the drift and the diffusion coefficients respectively for 
the optimum dynamics, and Wt is a Brownian motion which could be correlated or 
independent from of eq. ([T]). Based on these two OU evolutionary dynamics for 
the trait and its optimum, the model is called OUOU to reflect them. This model was 
established in [ 13 ] in order to study the adaptation between body size and tail length 
of woodcreepers. A special case of the OUOU model is the work in 0 which considers 
that the trait follows the OU dynamics of equation ([T]), and the optimum of the response 
trait, evolves via a Brownian motion, i.e. af = 0 and erf is time invariant. The OUBM 
model in turn was applied to study whether the sexual size dimorphism increases with 
body size when the female is the smaller sex in primates. 

The rate of evolution a\ in equation © measures the changing speed of the trait 
during its evolutionary process. A trait yt described in equation ([ 1 ]) may have a large 
variation when a species has evolved through the entire evolutionary history. Therefore 
this behavior cannot be captured by considering a constant rate of evolution as it has 
been presumed in the current existing literature. Indeed, there are many traits from a 
groim of related species with wide range of the evolutionary rate. For instance, Yopak et 
al. [3^ studied the variation in the brain organization of sharks. The widespread variation 
indicates the significant evolutionary diversity in their brain size and body mass. In such 
case, the variation should be captured in a stochastic way. Adopting these considerations, 
we develop and study in this paper the so-called OUBMBM and OUOUBM models, 
respectively, by treating af as a Brownian motion with a constant variance coefficient r, 

da\ = TdWr, (3) 

where the acronym OUBMBM (respectively OUOUBM) reflects the OU dynamics for 
the trait, a Brownian motion (respectively OU) for the optimum and a Brownian motion 
for the rate of evolution. Table 1 summarizes the models of trait evolution building on 
the general OU process that described by the SDE system of equation ([4]). 
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Table 1. Models for trait evolution developed under the OU process. Parameters in the table were either assumed 
as of constant values (0 is included if the model does not use the parameter), piecewise constants value or random 
variable. For instance, the OU model In [l3l described the trait evolution using the typical Orsntein-Uhlenbeck 
process while the OUmva model described the trait evolution under the generalized Hansen model with multiple 
optimum( 6 .y), multiple rate of evolution((T.y) and multiple constraining forces(a-,,), 7 = 1, 2, ■ ■ ■ , m. The BMv model 
assumed trait evolved with multiple rate of evolution under the Brownian motion [j^l . 


Section 2 presents preliminary results for our models of adaptation evolution and gives 
a precise mathematical formulation. Section 3 considers the novel models and hts them 
to large datasets in literature and compare them using the coefficient of determination 
(r-squared value), the Akaike information criterion(AIC) and an assessment of the bias 
of parameters. Last, Section 4 offers a discussion along with concluding remarks of this 
paper. 


2. Modeling adaptive evolution with random rate of evolution 

Adopting equations m-m, the trait evolution is organized in a vector form, Zj = 
{yt,0^, a^y, which satishes the following stochastics differential equation (SDE), 

dZt = AZtdt + DidWt, (4) 

0 \ 

where the upper triangle matrix At = I 0 —of 0 I consists of the parameters of 

\ 0 0 0 / 

selection strength; and Ct = DtD^ = diag((crt^)^, (cjf )^, r^) represents the associated 
covariance matrix and Wt = {W^, Wf, Wy)'^ is the vector of the associated independent 
Brownian motions. Note that the different noises are assumed mutually independent. We 
further assume that the force parameters are time invariant (i.e. = A is a constant 

matrix) and the rate of evolution for the optimum in equation ([ 2 ]) is a constant (i.e. 
cjf = ae). The SDE system described by Eq. (j4|) has a unique solution 

Zt = e-^'Zo + [' (5) 

Jo 

where Zq = (yo, Oq, is the initial condition for Z^ at t = 0. The expected value of the 
solution Zt, IE[Zf] = Zoe“^* and the second moment of the random vector Zt, denoted 
by Pi = E[ZiZ^], can uniquely be determined by solving the system of an ordinary 
differential equation 



APt + PiA^ + EQ. 


( 6 ) 
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Figure 1. An evolutionary tree for two species. The plot in the right panel shows the hypothetical evolutionary 
relationship between species i and species j . The plot in left panel represents the simulated trait evolution for the 
two species. Starting from time t = 0 (with trait value yo) to time t = ta (with trait value ya), the two species 
share a common ancestor. Then they diverge and evolve independently for t > ta into two species f , j at time 
t = -\- ti and t = ta -\- tj with trait value yi and yj , respectively. 


Next, we generalize the above model described by the SDE system in equation ([5]) for 
a group of n interacting species which share an evolutionary history described by a 
phylogenetic tree Let us consider that the current observed trait of the zth species, 
Ui^t, is a solution of the SDE given in equation © where 0 < t < T and T is the 
evolutionary time from the root of the tree to present time. Moreover, it is assumed 
that any pair of species, (i,j) share a most recent common ancestor at time instant 
t = ta where ta € [0, T] represents the evolutionary time from the root to the most 
recent common ancestor of the two species. To derive the joint distribution of the pair of 
random variables <t<T,we need to incorporate the corresponding shape 

and length (evolutionary time) from the phylogenetic tree into the models described 
in equation ([T]). Eigure [D displays a cartoon of a hypothetical evolutionary relationship 
between species i and species j where the affinity between the two species scaled in time 
unit can be represented by the following matrix G as shown in (l^ 


G = 


species i 
species j 


species i 
ta T ti 

ta 


species J 

ta 

ta T tj 


( 7 ) 


Let us dehne that E[yi\ya] is the expected trait value of species i conditioned on its 
ancestral trait value ya, where ya is the trait value of the most recent common ances¬ 
tor for species i and species j. The associated covariance of a pair of traits {yi,yj) for 
species i, j that diverged at time ta and evolved independently thereafter is given by 
Cov[yi,yj] = Cov[E[yi\ya], E[yi\ya]] [lli. The study in a considers a Brownian trait 
evolution which yields that the covariance of trait between species i and species j is 
proportional to the evolutionary time, that is Cov[yi,yj] = ayta- When an Ornstein- 


Uhlenbeck process is considered for a single trait evolution, the manuscript [l]| estab¬ 
lishes that the associated covariance, Cov[yi,yj] = j(2ay') when the initial condi¬ 

tion yo is assumed to be a random variable. Eurthermore, if conditioned yo on the root, 
then the associated covariance between two species trait evolved under the OU process 
is Cov[yi,yj\ = cjye“*o(i _ e“^““*“)/(2a) where the term tij denoted the evolutionary 
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distance since the two species diverged and evolved independently (i.e. tij = U + tj) (l^ . 
When the optimum of the response trait evolves randomly, Hansen et al. [l3l | and Jh- 
wueng and Maroulas 0 demonstrated that the variance-covariance structure between 
species i and species j under the OUBM model and OUOU model, respectively, can be 
derived under the following setting, 

Cov[yi,yj] = cl^Var[ya\ + cj^Var[ea\ -b 2cy_^ce^Cov[ya,0a], (8) 

where Var[ya],Var[6a] and Cov[ya,0a] can be derived from by solving equation 
(j6|) with initial condition at t = ta- In particular, we have Cy^ = and 

ce = for OUBMBM model and c„ = and cg = ———— 

for OUOUBM model. We provide the derivation of Cy^ and cg^ for the 
OUOUBM model in the Supplementary material. The Cy^ and cg^ for OUBMBM model 
can be derived in a similar manner. 

Once the covariance between two species is determined, the next step is to develop 
the statistical model for regression analysis for studying the adaptive relationship among 
traits. First, we quantify the evolutionary relationship between the predictor, xj, and 
the response trait, yt- Typically, the optimum of the response trait, has a linear 
relationship with the predictor xt (i.e. 9t = bo + biXt). Note that a general functional 
relationship between 6t and xt (i.e. dt = f{xt)) can indeed be the case and under this 
circumstance a linear approximation may be considered therein for parameter estimation. 
However, since the trait is typically transformed into the log scale before proceeding with 
the data analysis and the log transform converts nonlinear relationship into linear one, 
we use a linear relation in this work. 

Let {xi,yi),i = 1, 2, • • • , n be the pair of the trait values of species i observed at the 
tip of the phytogeny. The joint distribution of the entire evolutionary history for these 
two traits can be modeled by a bivariate random variable {xi^t,yi,t),0 < t < T. Since 
xt is assumed to have a linear relationship with the optimum of yt, the covariance for 
can be determined by dOt = bidxt- The rate of evolution for the optimum, 
ag, is identified once bi and cr^ are known( i.e. = bia'^). The predicted evolutionary 
regression of y on x can be derived as F^[y|x] = k + p{t)bix, where A; is a constant and 
p(t) = Cov[yt,9t]/Var[et], see e.g. 0. 

If the response Y = {yi,y 2 , ■ ■ ■ ,yn) and the predictor X = (xi,X 2 , • • • ,Xn)^ are ob¬ 
served at the tips of the phytogeny with the assumption that the relationship of the 
primary optimum 9t to the predictor variable is a simple linear regression, the evolution¬ 
ary regression curve for {X,Y) has the form 

y = X/3 + r (9) 

where X = [l,/9(t)Xi, • • • ,p{t)Xq] is the designed matrix of size 6 x y, 1 is the vector 
of ones, b = {bo,bi, ■ ■ ■ is a y dimensional vector of regression parameters, r is 

the residual vector following a normal distribution with zero mean vector, and with the 
residual covariance matrix V given below 

Yij = Cov[ri,rj] = Cov[yi - E[yi\9i],yj - E[yj\9j]] (10) 

where E[yt\9t] = /3o(t)-|-/3i(t)0t is the regression of the trait on the optimum (see Lemma 1 
in 0 for the derivation under the OUOU model, the optimal regression for other models 
described here can be derived in a similar manner). Note that equation (llOji involved four 
terms, Cov[yi,yj],Cov[yi, E[yj\9j]],Cov[yj, E[yi\9i]] and Cov[E[yi\9i], E[yj\9j]] which are 
computed once Sj was determined by solving equation ([6]) with initial condition Zq. 
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The statistical model derived from the OU dynamic of evolution is a multivariate 
normal distribution(i.e. Y ~ MVN(£i[X6], V)). The log likelihood for the regression 
analysis is 


logL(6,V|X,y,T) 


log 


^ ^-1{Y-Xbyv-^{Y-Xb) 

y/ (27r)" det(V) 


( 11 ) 


where T is a rooted phylogenetic tree with known topology and branch lengths and can 
be transformed to the matrix G directly for further use of constructing the matrix V. 
The latter is computed for the different models using the open source SAGE 28|. For 
parameter estimation, we use a similar algorithm as in as follows. 

Given the trait data Y and phylogenetic tree, the algorithm starts the search with an 
ordinary lease square estimate ho = (X^X)“^X^y. As T is given, the variance covariance 
structure for the residual V is calculated using equation (fTOt) . The likelihood in equation 
m is then optimized and the MLEs are recorded. The regression estimates b is then 
updated by 6 = (X V" X)"iX V" Y. This procedure is repeated until a distance be¬ 
tween the updated regression estimate and previous estimate is within an upper bound. 
We use the R packages optim 2^ for optimizing the likelihood. To address the potential 
sensitivity of the algorithm in the starting point, we use an alternative search for the 
MLEs where at most five different starting points are randomly selected in the domain of 
the parameter space and within each search, at most five attempts are set to access the 
convergence of estimation. If the convergence is not detected after the maximum itera¬ 
tion is reached, a different starting point will be either chosen using the current estimate 
or randomly selected in the domain. We keep updating the searches whenever a new 
improvement of the likelihood is observed. The search stops and save the estimates when 
either three improvement of likelihood are found or the maximum number of searches 
are reached. 


3. Empirical study and simulations 
3.1 Motivation 

Before any actual data set is analyzed using the models developed herein 
and their comparison with the current literature, we start by simulating 
the evolutionary path for the models. The paths were generated using pa¬ 
rameter values {ay,ay,ag,bo,bi) = (0.05,0.01,0.32,0.50,0.32) for the OUBM 

model; {ay,a 0 ,bo,bi,T) = (0.05,0.32,0.50,0.32,0.01) for the OUBMBM model, 

{ay,aQ,ay,a 0 ,bo,bi) = (0.01,0.01,0.20,0.45,0.05,0.30) for the OUOU model; and 
{oy, a 0 , ay, 00 , bo, 5i, r) = (0.01,0.01,0.20,0.45,0.05,0.30,0.01) for the OUOUBM model. 
Figure [3T] shows simulation paths under the OUBM/OUBMBM models (the left panel) 
and the OUOU/OUOUBM models(the right panel). We calculate the standard devia¬ 
tion (sd) of the trait considering the evolution paths from time t = 0 to t = 10000 in 
this simulation for each model. In this simulation, the OUBM model has sd of value 
72.64 while the sd is 107.96 for OUBMBM; The OUOU model has sd 1.61 while the 
sd is 8.00 for OUOUBM. The models OUBMBM and OUOUBM with the random rate 
of evolution account more variation than the special case models OUBM and OUOU, 
respectively which implies that evolutionary rate with high volatility should adopt a 
dynamic rate of evolution that is either an OUOUBM or OUBMBM model. 
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Figure 2. Evolutionary paths for the trait models. The path for the predictor trait xt(in black color) was first 
simulated under the Brownian motion or OU process; and the optimum paths Ot (in red color) was then simulated 
using the linear relationship Ot = bo + bixt; finally, the path for the response traits were simulated. In the left panel, 
yit represented path(in blue) simulated under the OUBM model, and y 2 t represented path(in purple) simulated 
under the OUBMBM; and in the right panel yn represented path (in blue) simulated under the OUOU model, 
and y 2 t represented the path(in purple) under the OUOUBM model. 


Table 2. 


Model Regression Line value 

OUOU y = 0.92 + 0.41a: 70^ 

OUBMBM y = 0.81 + 0.40x 72% 

OUOUBM y = 0.97 + 0.37a: 75% 

Evolutionary regression lines 


3.2 Shark dataset 


We first apply the proposed models herein that is the OUBMBM and OUOUBM as well 
as the OUOU model, which was suggested in 17| as the general version of 1^, in order 
to analyze the shark {chondrichthyans) dataset in [s^. A correlation study is conducted 
through the ordinary regression analysis and independent contrast method in and a 
signihcant relationship is found between the body size and brain mass. Figure [3] shows 
the evolutionary relationship represented by a rooted phylogenetic tree among the 42 
species of study interest. Due to high volatility in brain size (response variable) and 
the body weight (predictor variable), both datasets were log-transformed prior to data 
analysis. Figure [4] and Table shows the regression results. 

Comparing to the ordinary regression analysis(OLS) where y = 0.93 -|- 0.54x with 
= 74%, we find that the regression slopes for the four models were slightly shallower 
than the slope using the OLS approach in this dataset. On the other hand, the regression 
slopes for the models deviate from 0, therefore our models support the conclusion that the 
relative brain development reflects the dimensionality of the environment prey caption 
in addition to phylogeny in adaptation aspect as well as using the OLS and independent 
contrast method |35l | . In addition, we find that for this wide spread data the general 
OUOUBM model (r-squared value 75%) and the OUBMBM model (r-squared value 
72%) provides better fit than the OUOU model (r-squared value 70%). These results are 
summarized in Table [321 
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-Pristiophorus cirratus 
-Squalus sp 
-Squalus megalops 
-Squalus acanthias 
-Dalatias licha 
-Notorhynchus cepedianus 
-Carcharhinus brachyurus 
-Carcharhinus plumbeus 
-Carcharhinus leucas 
-Carcharhinus amblyrhynchos 
-Carcharhinus falciformes 
-Carcharhinus melanopterus 
-Thaenodon obesus 
-Negaprion acutidens 
-Prionace glauca 


-Sphyrna zygaena 
-Sphyrna mokarran 
-Sphyrna lewini 


-Galeocerdo cuvier 
-Hemigaleus microsstoma 
-Galeorhinus galeus 
-Mustelus lenticulatus 
-Mustelus canis 
-Mustelus antarcticus 
-Cephaloscyllium laticeps 
-Cephaloscyllium isabellum 
-Asymbolus rubiginosus 
-Asymbolus analis 


-Galeus boardmani 
-Isurus oxyrinchus 
-Carcharodon carcharias 
-Alopias vulpinus 
-Alopias superciliosus 
-Pseudocarcharias kamoharai 
-Carcharias taurus 
-Nebhus ferrugineus 
-Hemiscyllium ocellatum 
-Chiloscvilium punctatum 
-Brachaelurus waddi 
-Orectolobus maculatus 
-Orectolobus ornatus 
-Callorhinchus milii 


Figure 3. Evolutionary relationship of the sharks replotted from [^. 


3.3 Models comparison 


3.3.1 Accessing the statistical fit of the models via r-squared values 

^ Jhis Section, we consider 225 bivariate datasets appeared in the existing literature 
3" 21, 27, 3C, H, 341 and we compare the performance and fitting of the OUBMBM 


and OUOUBM models with the general OUOU model [17|]. First, we summarize the com¬ 
parison of the fit of these models using the coefficient of determination (r-squared values). 
The output of the analysis is given in Figure [S] where each point in the plot represents the 
r-squared value for the models of our interest. We compare the OUOU model with the 
OUBMBM model and summarize the result on the left panel. The comparison between 
OUOU model and OUOUBM model is given on the right panel. In each plot, the 1:1 line 
shows the equivalent fit of the both models; points below the 1:1 line indicate a better 
fit of the model whose r-squared is in the horizontal axis (OUOU model) and a point 
above the 1:1 line yields a better fit for the model that has r-squared in the vertical 
axis (the OUBMBM model and OUOUBM model). Overall both figures indicate that the 
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log body mass (g) 


Figure 4. Evolutionary regression curves. 


fit assessed by r-squared is consistent between the OUOU model and the new models 
(OUBMBM or OUOUBM) in most datasets. Most points are closed the diagonal line 
which indicates that when a high/low r—squared is observed by the OUOU model, the 
new model can identify a similar r—squared as well. Moreover, the OUBMBM model fits 
better than the OUOU model in 49.21 %. This shows that there might be no significant 
difference for two models as they have a samilar ability to detect the fit. On the other 
hand, the OUOUBM models proposed herein fits better than the existing OUOU model 
more frequently, in fact 57.40 % while the OUOU model hts better than the OUOUBM 
model for 42.60 % which accounts for a 14.80 % difference. 

3.3.2 Accessing the statistical fit of the models via the Akaike Information Criterion 

Next we compare the models employing their corresponding AICc values where 
AICc = —21ogL + 2nk/{n — k + 1), logL is the associated log-likelihood value, k is 
the number of model parameters and n is the number of extant species. We display the 
result in Figure [6] where the relative support for the models are reported by their corre- 
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Figure 5. Comparison of the models using r-squared values. 


spending Akaike weight u = exp(—O.SAA/Cc) where AAICc is the difference between 
the AICc value of the model and the minimal AICc value among the model set. Each 
horizontal line in Figure [6] represents the scaled Akaike weights of the models. 

For most datasets, we found that OUBMBM model accounts for more Akaike weight 
and happens to be the AICc selected models when compared to the OUOU model and the 
OUOUBM model. Since the OUOU model and the OUOUBM model only differ at most in 
one parameter from the OUBMBM model, this result can be resulted from the likelihood 
values where the OUBMBM contributes to lower likelihood than the OUOU/OUOUBM. 

3.3.3 Bias of parameters for the models 

We further access the bias of parameters for the models of interest. We also include the 
OUBM model in this section. We considered to generate data with four different sample 
sizes of n = 16,32,64, and 128 using the parameters values of {ay, ay) = (0.05,0.10) for 
OUBM, {ay,ao,ay) = (0.05,0.12,0.10) for OUOU, {ay,T) = (0.05,0.30) for OUBMBM 
and {ay,a 0 ,T) = (0.05,0.12,0.30) for OUOUBM. The true regression estimate {bo,bi) 
are set to (1.20,0.72), and 50 replicates are sampled from this set up with considering 
of using four different types of tree phylogeny (1) a star tree where all species are com¬ 
pletely unrelated, (2) a completely balanced tree, (3) a completely pectinate tree, and 
(4) a random tree generate under the birth-death process. Hence 4x4x4x 50 = 800 
replicates are generated and the associated parameters are estimated through the MLE 
approach described in Section 2. For each parameter of interest, we report the boxplots 
of parameter estimates ay,ag,ay, and r under each models and different sample sizes 
(we combine the tree topology effect in this study). The results are shown in Figure Sl- 
S4 (see supplemental material section). Overall, we find that the accuracy for estimating 
the parameters of interests is improved with increasing the sample size. The interquartile 
range(IQR) shrinks with larger sample size increases for each model. 

Figure SI and Figure S2 shows that estimation for Oy and ag, respectively, cannot 
achieve satisfactory accuracy because there are many outliers configured in the corre¬ 
sponding boxplots. Focusing on Figure SI, the parameter ay is in general poorly es¬ 
timated under the OUBM and the OUBMBM model. This situation also occurs for 
estimating ag as in Figure S2 the parameter estimates become less spread as the sample 
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ouou oubmbm ououbm 



AlCc weight 


Figure 6. The relative support measured by AICc weight among the three models in this study (OUOU, 
OUBMBM, and OUOUBM) across 225 datasets where each horizontal line represents a weight of the models. 
The OUBMBM model accounts for more support for a majority of the datasets. 


size increases. Figure S3 and Figure S4 show that parameter estimation for ay and r, 
respectively, are in general good. In Figure S3, OUOU model accounts for larger varia¬ 
tion in estimating Gy than the OUBM model. Furthermore, Gy for the OUOU model can 
be estimated more accurately as the sample size increases. In Figure S4, the OUOUBM 
model provides wider IQR than the OUBMBM model for estimating the parameter r. 
We leave a discussion for this in the next session. 


4. Conclusion 

In this paper, we developed two models for the adaptive trait evolution and evalu¬ 
ated their performance by analyzing many empirical datasets. We found that our model 
OUBMBM/OUOUBM fit better for more datasets than their submodel model OUOU 
when evaluated the performance of fit by the r-squared values. On the other side, under 
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the likelihood based model selection criterion, we found that due to the likelihood, the 
OUBMBM model became AICc selected models for most empirical datasets than the 
OUOU and OUOUBM models. 

In bias study of parameter, we found that the parameters can not be always estimated 
accurately under the MLE analysis for many datasets. Ho and Ane 0 and Ane [1 
described the limitation of the parameter estimation of the BM model described in 91] 
and the limitation of the OU model described in 0 for trait evolution. They pointed 
out that since some parameters may not identifiable, the maximum likelihood estimators 
for trait models could fail to be consistent estimators where the convergence to the true 
parameter cannot be guaranteed. For the OU model. Ho and Ane 0 proved that the 
selection optimum cannot be estimated consistently as the tree grows indefinitely. In 
our framework, we might encounter the same problem when estimating the parameters 
for the models of adaptive evolution developed here. To deal with this situation, we 
suggest a future study for the models of adaptive evolution under a Bayesian paradigm 

for works 
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to reduce the estimation difficulty for non-identifiable parameters (see 
of comparative methods using a Bayesian approach). 

The simulation paths in Figure 13.11 indicated the higher variation of the OUBMBM 
model with respect to the OUBM model, and the higher variation of the OUOUBM 
model with respect to the OUOU model. It is expected that the OUBMBM/OUOUBM 
model would provide a better fit than the OUOU model for data with wider variation. 
In the shark dataset, we demonstrate that the both OUBMBM and OUOUBM models 
provide a better fit than the OUOU model. 

The SAGE source code for solving the moments of and the covariance between the 
residual V, the R source codes for parameter estimation and simulations and the datasets 
used in this work can be accessed directly at www.tonyjhwueng.info/OUOUBMsim. 


4.1 Acknowledgements 

4.2 Funding 

4.3 Supplemental material 

4-3.1 Variance-covariance structure for the OUOUBM model 

!—Oiy ay 0\ 

For OUOUBM model, A = I 0 —ag 0 1 . Solving equation ([5]) amounts to calculate 

\ 0 0 0 / 

the exponential A. Note that A has three distinct eigenvalues {—ay, —ag,0} with eigen¬ 
vectors set Q = {(1,0,0)', 1,0)', (0,0,1)'}. This yields that the term can be 

computed straightforwardly by e~^^ = QAQ~^ where A = diag{e“"“^ 1} is the 

diagonal matrix and Q = {vi,V 2 ,V 3 } is the matrix of eigenvectors. Since the white noises 
Wt have expectation zero, the expected value of can be determined given the initial 
condition Zq = {yo,0o, cro)'. In particular, the expecated value of yt conditioned on the 
initial value y(0) = yo has the form 


E[yt\yo] = Cy,yo + cg^Oo 


( 12 ) 


*andce, = - e “A). 

Then we can apply equation m to yt conditioned on any ancetral value, that is 


where Cy^ = e " 


E[yt\ya] = Cyj?/a + where ya and da are the ancetral value at t In [l]|, the 

associated covariance of a pair of traits {yi,yj) for species i, j that diverged at time 
ta and evolved independently thereafter is given by Cov[yi,yj] = Cov[E[yi\ya], E[yi\ya]] 
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Figure SI. Boxplot for cty. 

where E[yi\ya] = + p, 9a and Uj is the evolutionary 

distance since species i and i diverged. 

To complete the calculation of Cov\yi,yj\, the next step is to calculate the terms 
Cov{ya,ya),Cov{ya,6a) and Cov{9a,9a). These term can be determined by solving the 
ordinary differential equations (ODEs) in equation ([6]). Notice that although there are 
3x3 = 9 ODEs in Eq. ([6]), due to symmetry of Pt it suffices to solve six equations ( 
including the three equations for the second moments of yt, 9^, at) and the expected value 
of yt9t,yt(^t and 9tat). Since some variables are embbeded in the equations, we cannot 
solve the six ODEs simultaneously. Fortunately, we can determine a solution recursively 
once upon one ODE is completely solved and it is then used for solving another ODE. 
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